This invention is in the field of satellite navigation. Embodiments of this invention are directed to Global Navigation Satellite System (GNSS) receivers, such as of the Global Positioning System (GPS) type, and to improved estimation of position and velocity obtained by such receivers comprehending measurement uncertainties.
In recent years, satellite navigation equipment has become widespread. Sophisticated navigation systems for large-scale vehicles such as aircraft and ships, and positioning systems for map and survey functions, were the first implementations of GNSS technology. In recent years, however, GPS receivers have been implemented in relatively modest applications, including many consumer automobiles, add-on navigation systems for automobiles, and handheld navigation systems for hikers. Indeed, GPS systems are now included in golf carts to provide golfers with accurate distance information, and are now commonplace in mobile telephone handsets.
GNSS technology, currently implemented by way of the GPS system, has many advantages over previous navigation systems, particularly land-based systems such as LOng-RAnge Navigation (LORAN) systems. GPS signals can be used at all times of the day and night, and are not greatly affected by weather or atmospheric conditions. Line-of-sight transmission and worldwide coverage are facilitated by the number and positions of transmitting satellites in orbit. The GPS system also provides both position and velocity information to the user. Varying degrees of position and velocity accuracy are available, enabling a large number of users of inexpensive receivers, while still allowing sophisticated users to obtain a higher level of accuracy.
In a general sense, by way of background, GPS navigation is based on triangulation of the receiver location using signals received from multiple satellites of known location, in which the received signals are effectively time-stamped with the time of transmission. More specifically, GPS satellites periodically transmit a pseudo-random number (PRN) sequence via a spread-spectrum signal, in which the transmitted value of the frame boundaries of the data modulating the PRN sequence, as well as the PRN sequence itself, has a deterministic relationship to a known time (i.e., the beginning of the GPS “week”); as such, the received data from the satellites effectively include a time stamp. From the PRN value, GPS receivers determine the signal propagation time from the satellite to the receiver and, multiplying that propagation time by the speed of light, determine a measured distance (the “pseudorange”) of the receiver from each of the multiple satellites. Typically, the pseudoranges from four or more satellites of known position, in an earth-centered earth-fixed coordinate system, are then triangulated by solving a system of position equations, to determine the position of the receiver in that coordinate system.
FIG. 1 illustrates, in an idealized 2-D representation, the fundamental concept of GPS triangulation. In this illustration, each of satellites SAT1, SAT2, SAT3 have transmitted a time-stamped signal that is received by a synchronized receiver at an unknown location. By subtracting the time of transmission indicated in the signal time stamp from the time of receipt, and multiplying that difference by the speed of light, the receiver can estimate its distance r1, r2, r3 from respective satellites SAT1, SAT2, SAT3 of known position within the coordinate system being used. Computational circuitry within the receiver can then identify point RLOC in that coordinate system at which all of the estimated distances r1, r2, r3 coincide, to a best approximation. In the 3-D sense, each radial distance r1, r2, r3 from satellites SAT1, SAT2, SAT3 will of course define a sphere, requiring either four or more satellites or extrinsic information (e.g., knowing that location RLOC is on the surface of the earth) to resolve a unique receiver location RLOC.
As known in the art and as mentioned above, the geometric “range”, or distance, r from a GPS satellite to a GPS receiver corresponds to the propagation time between the transmission time Ts at the satellite and the receipt time Tu at the user, multiplied by the speed of light c:r=c(Tu−Ts)=cΔt The “pseudorange” that can be determined by the GPS receiver necessarily involves consideration of the time offset tu between the time at the receiver and the true reference time (“system time”), and the time offset δt between the time at the satellite and system time. As such, the pseudorange ρ considering these time offsets can be expressed as:
                    ρ        =                ⁢                  c          ⁡                      [                                          (                                                      T                    u                                    +                                      t                    u                                                  )                            -                              (                                                      T                    s                                    +                                      δ                    ⁢                                                                                  ⁢                    t                                                  )                                      ]                                                  =                ⁢                              c            ⁡                          (                                                T                  u                                -                                  T                  s                                            )                                +                      c            ⁡                          (                                                t                  u                                -                                  δ                  ⁢                                                                          ⁢                  t                                            )                                                              =                ⁢                  r          +                      c            ⁡                          (                                                t                  u                                -                                  δ                  ⁢                                                                          ⁢                  t                                            )                                          Or, in coordinate space, defining s as the vector from the origin (i.e., center of the earth, in an earth-centered coordinate system) to the satellite, u as the vector from the origin to the receiver location, and r as the vector from the receiver location to the satellite corresponding to the vector difference s−u:ρ−c(tu+δt)=∥s−u∥Pseudorange ρ and the time offset values are measurable or otherwise knowable, thus allowing solution of the actual range. And considering that the GPS network itself now determines corrections for the satellite clock offset δt from system time, and transmits those corrections to the satellites for rebroadcast to receivers, satellite clock offset δt is effectively compensated for at the receiver, such that:ρ−ctu=∥s−u∥
Conventional GPS systems determine the user position (i.e., receiver location RLOC in FIG. 1) based on pseudorange measurements to four or more satellites, giving rise to a system of equations (one for each satellite):ρj=∥sj−u∥+ctu where j is the satellite index. For a system of four satellites (j=1, 2, 3, 4) in which satellite j is at position (xj, yj, zj) in the coordinate system, the pseudoranges ρj can be expressed into a system of equations in which the unknowns include the user (receiver) position (xu, yu, zu) and the user time offset tu:ρ1=√{square root over ((x1−xu)2+(y1−yu)2+(z1−zu)2)}{square root over ((x1−xu)2+(y1−yu)2+(z1−zu)2)}{square root over ((x1−xu)2+(y1−yu)2+(z1−zu)2)}+ctu ρ2=√{square root over ((x2−xu)2+(y2−yu)2+(z2−zu)2)}{square root over ((x2−xu)2+(y2−yu)2+(z2−zu)2)}{square root over ((x2−xu)2+(y2−yu)2+(z2−zu)2)}+ctu ρ3=√{square root over ((x3−xu)2+(y3−yu)2+(z3−zu)2)}{square root over ((x3−xu)2+(y3−yu)2+(z3−zu)2)}{square root over ((x3−xu)2+(y3−yu)2+(z3−zu)2)}+ctu ρ4=√{square root over ((x4−xu)2+(y4−yu)2+(z4−zu)2)}{square root over ((x4−xu)2+(y4−yu)2+(z4−zu)2)}{square root over ((x4−xu)2+(y4−yu)2+(z4−zu)2)}+ctu Theoretically, of course, this system of four equations and four unknowns has a unique solution. In practice, however, each pseudorange measurement includes some amount of noise or error. In the presence of such noise, according to conventional techniques, this system of nonlinear equations can be solved, to derive user position (xu, yu, zu) and user time offset tu, by way of iterative techniques such as least squares minimization applied to a linearization of these non-linear equations, or through the use of Kalman filtering over a time sequence of the pseudorange measurements.
It is useful, for purposes of comprehension, to describe the operation of an example of the linearization approach to the solution of this system of equations. As known in the art, linearization can be carried out with reasonable accuracy from an approximate position ({circumflex over (x)}u, ŷu, {circumflex over (z)}u), for example as may be determined by cellular positioning or other known methods. In the absence of assistance (i.e., via “assisted GPS” as known in the art), the initial approximate position can be taken as the center of the earth; as such, the approximate position need not be very accurate. The offset of the true user position (xu, yu, zu) from this approximate position can then be denoted by a displacement (Δxu, Δyu, Δzu). Taylor expansion of each of the set of pseudorange ρj equations allows this displacement or position offset (Δxu, Δyu, Δzu) to be expressed as linear functions of the known coordinates and the pseudorange measurements. With a pseudorange ρj represented by:
                                          ρ            j                    =                    ⁢                                                                                          (                                                                  x                        j                                            -                                              x                        u                                                              )                                    2                                +                                                      (                                                                  y                        j                                            -                                              y                        u                                                              )                                    2                                +                                                      (                                                                  z                        j                                            -                                              z                        u                                                              )                                    2                                                      +                          ct              u                                                                    =                    ⁢                      f            ⁡                          (                                                x                  u                                ,                                  y                  u                                ,                                  z                  u                                ,                                  t                  u                                            )                                            ⁢        one can calculate an approximate pseudorange {circumflex over (ρ)}u using the approximate position ({circumflex over (x)}u, ŷu, {circumflex over (z)}u) and a time bias estimate {circumflex over (t)}u:
                                          ρ            ^                    j                =                ⁢                                                                              (                                                            x                      j                                        -                                                                  x                        ^                                            u                                                        )                                2                            +                                                (                                                            y                      j                                        -                                                                  y                        ^                                            u                                                        )                                2                            +                                                (                                                            z                      j                                        -                                                                  z                        ^                                            u                                                        )                                2                                              +                      c            ⁢                                          t                ^                            u                                                              =                ⁢                  f          ⁡                      (                                                            x                  ^                                u                            ,                                                y                  ^                                u                            ,                                                z                  ^                                u                            ,                                                t                  ^                                u                                      )                              Each component of the true user position (xu, yu, zu) and receiver clock offset tu can be expressed as the sum of an approximate component and an incremental component:xu={circumflex over (x)}u+Δxu yu=ŷu+Δyu zu={circumflex over (z)}u+Δzu tu={circumflex over (t)}u+Δtu in which case:ƒ(xu,yu,zu,tu)=ƒ({circumflex over (x)}u+Δxu,ŷu+Δyu,{circumflex over (z)}u+Δzu,{circumflex over (t)}u+Δtu)This function ƒ(xu, yu, zu, tu) can be expanded about the approximate position and time bias estimate ({circumflex over (x)}u, ŷu, {circumflex over (z)}u, {circumflex over (t)}u) using a Taylor series and truncating the partial derivative terms after the first-order expressions, so that the pseudoranges ρj can be expressed as:
      ρ    j    =                    ρ        ^            j        -                                        x            j                    -                                    x              ^                        u                                                r            ^                    j                    ⁢      Δ      ⁢                          ⁢              x        u              -                                        y            j                    -                                    y              ^                        u                                                r            ^                    j                    ⁢      Δ      ⁢                          ⁢              y        u              -                                        z            j                    -                                    z              ^                        u                                                r            ^                    j                    ⁢      Δ      ⁢                          ⁢              z        u              +          ct      u      where {circumflex over (r)}j is defined as:{circumflex over (r)}j=√{square root over ((xj−{circumflex over (x)}u)2+(yj−ŷu)2+(zj−{circumflex over (z)}u)2)}Introducing shorthand variables axj, ayj, azj as follows:
            a      xj        =                            x          j                -                              x            ^                    u                                      r          ^                j                        a      yj        =                            y          j                -                              y            ^                    u                                      r          ^                j                        a      zj        =                            z          j                -                              z            ^                    u                                      r          ^                j            and defining pseudorange difference Δρ=ρj−{circumflex over (ρ)}j, one can express the set of linear equations to be solved, based on range measurements to four satellites, as:Δρ1=ax1Δxu+ay1Δyu+az1Δzu−ctu Δρ2=ax2Δxu+ay2Δyu+az3Δzu−ctu Δρ3=ax3Δxu+ay3Δyu+az4Δzu−ctu Δρ4=ax4Δxu+ay4Δyu+az4Δzu−ctu By expressing the coefficients and variables in matrix form:
            Δ      ⁢                          ⁢      ρ        =          [                                                  Δ              ⁢                                                          ⁢                              ρ                1                                                                                        Δ              ⁢                                                          ⁢                              ρ                2                                                                                        Δ              ⁢                                                          ⁢                              ρ                                  3                  ⁢                                                                                                                                                                Δ              ⁢                                                          ⁢                              ρ                4                                                        ]        ;      H    =          [                                                  a                              x                ⁢                                                                  ⁢                1                                                                        a                              y                ⁢                                                                  ⁢                1                                                                        a                              z                ⁢                                                                  ⁢                1                                                          1                                                              a                              x                ⁢                                                                  ⁢                2                                                                        a                              y                ⁢                                                                  ⁢                2                                                                        a                              z                ⁢                                                                  ⁢                2                                                          1                                                              a                              x                ⁢                                                                  ⁢                3                                                                        a                              y                ⁢                                                                  ⁢                3                                                                        a                              z                ⁢                                                                  ⁢                3                                                          1                                                              a                              x                ⁢                                                                  ⁢                4                                                                        a                              y                ⁢                                                                  ⁢                4                                                                        a                              z                ⁢                                                                  ⁢                4                                                          1                              ]        ;            Δ      ⁢                          ⁢      x        =          [                                                  Δ              ⁢                                                          ⁢                              x                u                                                                                        Δ              ⁢                                                          ⁢                              y                u                                                                                        Δ              ⁢                                                          ⁢                              z                u                                                                                                        -                c                            ⁢                                                          ⁢              Δ              ⁢                                                          ⁢                              t                u                                                        ]      the system of equations can be expressed as:Δρ=HΔx and, in theory, can be solved for the vector of position and time displacement Δx from the estimated position and time offset via:Δx=H+Δp where H+ represents a pseudoinverse of linearized derivative matrix H, calculated according to the particular solution technique (e.g., least-squares, weighted least-squares, etc.).
In similar fashion, GPS navigation provides the capability of determining the three-dimensional user velocity {dot over (u)} as an approximation of the time rate of change of the calculated user position u over a time interval:
      u    .    =                    ⅆ        u                    ⅆ        t              =                            u          ⁡                      (                          t              2                        )                          -                  u          ⁡                      (                          t              1                        )                                                t          2                -                  t          1                    Doppler shifts in the frequencies of the received satellite signals are analyzed to solve the velocity problem, in many modern GPS receivers. In a linearization approach similar to that described above for determining the user position, a variable dj is used for the jth satellite signal:dj={dot over (x)}uaxj+{dot over (y)}uayj+żuazj−c{dot over (t)}u For signals received from four satellites, a system of four equations in the four unknowns (i.e., the derivative terms {dot over (x)}u, {dot over (y)}u, żu, {dot over (t)}u) can be expressed in matrix form:d=Hg where gT=[{dot over (x)}u {dot over (y)}u żu {dot over (t)}u], permitting solution of vector g from this system via:g=H+d 
As mentioned above, in practice, these calculations by way of which GPS systems can determine a user position and a user velocity are corrupted by noise in the signal. Referring to user position, the measured pseudoranges ρ amount to a sum of the computed pseudoranges h(x), based on the true user position and time offset vector x, with a composite measurement noise vector v:ρ=h(x)+v For the least-squares solution method, the linearized system of equations including the effects of noise is more accurately expressed as:Δρ=HΔx+v and, in theory, can be solved for the vector of position and receiver clock bias Δx from the estimated position and time offset via:Δ{circumflex over (x)}=H+ΔρwhereΔ{circumflex over (x)}=Δx+e=(x−{circumflex over (x)})+e {circumflex over (x)} being the estimated user position and receiver clock bias, and e representing the position-domain error caused by the measurement noise v in the pseudorange measurements, as applied to the derivative matrix H:e=H+v 
Measurement noise vector v of the signals as received from the satellites depends on such factors as pure measurement noise, errors included within the applicable model considered by the receiver, and any unmodeled effects. In conventional GPS systems, to the extent that random variables are included within vector v, those random variables are typically modeled as uncorrelated with one another, and as having equal variances over satellite signals.
It is contemplated, however, according to this invention, that the equal variance assumption is invalid in many practical situations. For example, a GPS receiver located within an “urban canyon” will receive GPS signals from some satellites in highly attenuated form, if not completely blocked by the surrounding buildings. In this case, some of the satellite signals received by that GPS receiver may be reflected signals, or may be received over multiple paths because of the reflections (i.e., multipath signal transmission). According to this invention, it is contemplated that these suboptimal received satellite signals will have a measurement noise vector v that is time-varying due to the changing channel environment as the user moves through the physical environment. These signals will have unequal variances among one another, and will generally exhibit time-varying non-zero biases. As such, it has been observed, in connection with this invention, that treatment of the random variables in the measurement noise vector v as uncorrelated and having equal variances will lead to inaccurate position and velocity determinations.